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ABSTRACT 


We  consider  plane -wave  motion  at  normal  incidence  in  a horizontally 
layered  system.  The  system  is  assumed  lossless,  and  only  the  compress- 
ional  waves  are  treated.  A procedure  is  introduced  for  determining  the 
reflection  coefficients  of  the  layered  system  when  the  observed  seismic 
data  may  contain  random  noise.  No  deconvolution  of  the  measured  seismic 
data  is  required  by  the  procedure  when  the  input  is  a narrow  wavelet. 
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1.  INTRODUCTION 


In  recent  years  much  attention  has  been  given  to  the  problem  of 
determining  reflection  coefficients  for  a layered  media  from  the  observed 
seismic  data  [1-4].  In  line  with  the  customary  assumptions  and  restrictions, 
we  also  limit  our  attention  to  a horizontally  stratified  nonabsorptive  earth 
with  vertically  traveling  plane  compressional  waves.  Such  a system  is 
completely  described  by  a set  of  reflection  coefficients  and  travel  times 
within  layers. 

A fundamental  procedure  described  in  detail  in  the  above  references 
for  deriving  values  of  the  reflection  coefficients  can  be  summarized  by 
the  following  assumptions  and  steps. 

Standard  Assumptions: 

(A  1 ) The  input  wavelet  is  assumed  known. 

(A2)  The  data  is  assumed  noise  free. 

(A3)  The  layered  system  is  assumed  to  have  uniform  travel 
times  between  layers  where  a number  of  the  layers  are 
hypothetical,  i.  e.  , fliey  may  not  correspond  to  an  actual 
interface  of  the  substructure  and  are  associated  with 
zero  reflection  and  unity  transmission  coefficients. 

Standard  Steps : 

(51)  The  observed  seismic  data  is  deconvolved  using  the  input 
waveform  to  produce  the  system  response  to  a unit  spike 
input. 

(52)  The  number  of  layers  is  chosen  high  enough  to  result  in 
travel  times  short  compared  with  the  inverse  of  the 
bandwidth  of  the  observed  seismic  data. 

(53)  The  deconvolved  data  is  sampled  with  sampling  interval 
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equal  to  the  chosen  one-way  travel  time  between  layers. 

(54)  The  system  structure  is  used  to  arrive  at  a set  of 
normal  equations  (linear  simultaneous  equations)  in 
terms  of  reflection  coefficients  and  the  discretized  and 
deconvolved  observed  data. 

(55)  The  normal  equations  of  the  preceding  step  have  the 
Toeplitz  structure  which  makes  it  possible  to  utilize 
the  very  efficient  Levinson  algorithm  to  recursively 
solve  for  the  reflection  coefficients. 

In  this  paper  the  method  of  solution  to  the  inverse  problem  stated 
above  is  fundamentally  modified  to  cope  with  the  existence  of  the  noise  in 
the  measurement  data,  often  without  need  for  any  deconvolution.  More 
specifically,  although  again  a uniform  layered  system  is  assumed,  the 
choice  of  number  of  layers  can  now  be  made  independent  of  the  sampling 
rate  requirement  of  the  data  (step  (S2)  above)  often  resulting  in  the  need 
for  far  fewer  layers.  No  deconvolution  is  necessary  (step  (SI))  for 
wavelets  of  duration  of  the  order  of  twice  the  layer  travel  times.  The 
exact  deconvolution  of  step  (SI)  is  either  not  possible  in  practice  or,  at  the 
least,  will  further  aggrevate  the  harmful  effects  of  the  noise  in  the  obser- 
vation [5].  Furthermore,  the  deconvolution  is  a time  consuming 
operation.  Finally,  the  procedure  is  very  simple  to  derive  and  does  not 
need  the  concepts  of  z-transforms , minimum  phase,  forward  and  backward 
polynomials,  spectral  factorization,  etc.  The  results  reduce  to  the 
existing  solution  of  the  inverse  problem  in  the  absence  of  noise  and  with 
a spike  input  signal  (wavelet)  [l]. 


Of  course  the  deconvolution  may  be  performed  in  discrete  time  using 
the  same  sampling  interval. 
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2..  STATEMENT  OF  THE  PROBLEM 


We  are  considering  a uniform  K layered  system  and  normal 

incident  compressional  waves.  Figure  1 represents  such  a system 

th 

where  d.(t)  is  the  down-going  wave  at  the  bottom  of  the  j layer  and 
u.(t)  is  the  up-going  wave  at  the  top  of  the  layer.  The  reflection,  down- 
ward transmission  and  up-ward  transmission  coefficients  associated 

with  the  interface  at  the  bottom  of  layer  are  denoted  r.,  t.  and 

J J 

I I 

t.  respectively  where  t.=  l+r.,  t.=  l-r..  The  one  way  travel  time 
J J J J J 

between  layers  is  denoted  by  t . 


1 "o**' 

* 

c 

o 

^Uj(t) 

^Uj(t) 

iV.'*' 

layer  0 (half  space) 


layer  1 


layer  j 


layer  K+1  (Basement) 
Figure  1.  K Layered  System 


- 3 - 


T 


The  input  to  the  system,  ^^(t),  is  assumed  known  (the  wavelet) 
and  the  output  may  be  either  Uj(t)  (in  the  marine  environment)  or 
Up(t).  The  measured  seismic  data,  y(t),  consists  of  the  output  and  an 
additive  noise  component  n(t).  The  source  of  this  noise  may  be  the 
instrument  measurement  noise,  the  uncertainty  in  the  knowledge  of  the 
input  wavelet  or  response  to  unwanted  inputs  (ambient  noise).  It  is  desired 
to  process  y(t)  , t SO  and  derive  values  for  the  reflection  coefficients 
Tj  , j = 0,  1,...,K  (r^  may  be  assumed  known  in  cases  such  as  the 
marine  environment). 


3.  STATE  EQUATIONS 


Using  the  notation  of  Figure  1,  for  a general  layer  we  have 


[6,7]. 

+ = ‘j  + r.  d.(t)  ( 1 ) 

Vi 

These  equations  are  valid  for  j=l,...,K-l.  They  should  be 
augmented  at  the  surface  with 

V‘'  = + 

<i|(t  + T)  = Ui(l)  + Iq  do(t)  (4) 

and  at  the  basement  with 

“k<'  = 'k  Vtl'*'  * ■'k  ' "k  “k"'  < 5 > 

“Ktl*''  = -^K-l  “km'*'  **K  V'>  = ‘k  V*'  < ^ ' 
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^ — ' — 

Equations  (3,4)  and  (5,6)  can  be  derived  from  (I)  and  (2)  (^letting 

j = 0,  1,,..,K]  by  noting  that  UQ(t)  is  taken  at  the  bottom  of  layer  0 

and  d (t)  at  the  top  of  layer -(K+1)  and  n , (t)  0.  The  term 

tv+i  K+1 

dj^^j(t)  represents  the  down-going  wave  leaving  the  last  interface  and 
is  not  reflected  by  any  other  interface;  hence  u (t)  - 0.  These 
equations,  called  causal  functional,  are  not  difference  equations  since 
t is  the  continuous  time  variable  [6] . 


4.  A GENERALIZED  ENERGY  TRANSFER  (KUNETZ)  RELZ-TICN 

Consider  e to  be  a non-negative  continuous  or  discrete  variable 
with  dimension  of  time.  Equations  (1)  and  (2)  [wher  e j 0,  1 , . . . , K] 
are  multiplied  by  u^(t+T  + e)  and  d^^^(t+T  + e)  respectively  resulting  in 


Uj(t+T)  u^(t  + T+ e)  = t^^  u^^^(t)  u^^^(t+ e)  + r^^  d^(t)  d^(t+ e) 

+ r . t'.[u  (t)  d.(t  + e)  + u.  ,(t+E)d.(t)] 

J J+1  J J+1'  J ^ 


( 7 ) 


dj^l(t  + T)  dj^l(t  + T+E)  = r^^  u^^j(t)  u^^j(t+ e)  + t^^  d.(t)  d^(t+ e) 

- r.  t.[u.^^(t)  d.(t+  E)  + u^^j(t+  e)  d.(t)] 


( 8 ) 


Multiplying  (7)  by  t^/t^  and  adding  the  resulting  expression  to  (8)  yields 


j(t  + t)  d^^j(t  + t + e)  + + ’’’)  + T + e) 

~ d^(t  + e)  + Uj^j(t)  Uj^j(t+e) 


( 9 ) 
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Lee  us  define  the  following  correlation-type  functions 


r 


D,(e)  = r d.(t)  d.(t+  e)  dt  (10) 

J J J J 


A 

u.(£)  = i u.(t)  u.(t+ e)  dt  (11) 

J J J J 


Integrating  both  sides  of  Eq.  (9)  from  to  + “ , and  using  Eqs.  (10)  and 
(11),  we  find  that 


D.  ,,(e)  - U.  ,,  (e)  = t /t:[D.(£)  - U.(e)l  (12) 

J+1  J+1  J J J J 


where  j = 0,  1 , 2 , . . . , K.  This  is  a generalization  of  the  well-known  [ 1 ] 
energy  transfer  (Kunetz)  relation.  Note  that  in  our  derivation,  input  ^^(t) 
is  not  assumed  to  be  an  impulse  and  the  seismic  data  is  not  discretized. 
An  abstract  generalization  of  (12)  is  given  in  Appendix  A. 

Iterating  (12),  starting  with  j - £ and  ending  with  j K,  we  obtain 


K 

n t. 

n t; 
, 1 


where  i can  take  on  the  values  of  0,  1 , . . . , or  K.  In  the  marine  case  this 
relationship  is  used  with  2=1.  In  the  non-marine  case  (Appendix  B), 
it  is  used  with  j2  = 0. 
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5.  APPLICATION  TO  MARINE  ENVIRONMENT 


In  this  section  we  will  direct  our  attention  to  the  marine  case 

and  shall  express  D ,{e)  in  terms  of  measured  signals.  To  do 

Kt  X 

this,  we  set  r^  = 1 and  we  see  from  (13)  that  we  must  express 
Dj(£)  - Uj(e)  in  terms  of  the  measured  signals.  The  first  layer 
can  be  depicted  as  in  Figure  2. 


todJt)J  u,_(t) 


-Uj(t) 


Figure  2.  First  layer  in  the  marine  case. 


Observe  that  (4)  becomes 


dj(t  + T)  = -Uj(t)  f 2dQ(t)  . 


( 14  ) 


From  (10),  (11),  and  (M),  v.  c can  evaluate  the  diffcTence  term  D^(0-U^('), 


Dj(£)-Uj(e)  ^ P(t)  = J i[2dQ(t) -u^(t)l[2dQ(t+ e) -Uj(t+ E)] 

.CD 

- Uj(t)  Uj(t  + e))  dt 


P(E)  = f 4d^(t)  dQ(t  + E)  dt  - J 2d^(t)  Uj(t+  *)  dt 


r 2u  J (t)  d^(t  + e)  dt  ; 


( 15  )•■■ 


( 16  ) 


* Because  of  the  range  of  integration  in  (10),  we  can  also  express  D;(i)^is 

D.(k)=  I d.(t+T)d.(t+T  + f)dt.  We  use  this  form  of  (10)  in  our  de velopm cnt 
j . on  j j 


of  Dj(t)  - Uj(r). 


1 


hence,  P(e)  can  be  evaluated  from  a knowledge  of  dQ(t)  and  u^(t) 
for  any  desired  e.  Observe,  also,  that  (13)  with  can  be  written 

in  terms  of  F(f),  using  (15),  as 


K 

n t. 

1 ‘ 


K , 
n t. 
1 ‘ 


P(E) 


( 17  ) 


We  should  point  out  at  this  stage  that  the  quantity  Uj(t)  needed 
in  (1())  is  only  available  through  the  observation 

y(t)  = Uj(t)  + n(t) 

where  n(t)  is  the  additive  noise.  Consequently,  P(e)  is  not  physically 
available:  however,  we  can  define  P(0  by  replacing  y(t)  for  Uj(t)  in 


(lb). 


+CO 


P(E)  ^ J 4 djj(t)  dQ(t+ e)  dt  - j 2 dQ(t)  y (t+ e)  dt 


( 18  ) 


+ CD 


2y  (t)  dQ(t  + E)  dt 


which  can  also  be  written  as 

P(£)  = P(e)  + N(e)  , 

where 

+0O 

N(e)  = -'  2 d^(t)  n(t  + E)  dt  - ■ 2 n(t)  d (t  f e ) dt 

*i  0 *1  V 


( 19 


( 20  ) 


The  statistics  of  noise  term  N(e)  can  be  determined  in  terms  of 
these  of  n(t).  Using  P(r)  in  place  of  P(e)  in  (17)  yields 


8 


+ = I 'i  V 


( 22  ) 


Note  that  the  coefficient  of  the  highest  term  of  the  left  hand  side  is  unity 

and  that  of  the  lowest  term  is  r r The  precise  form  of  the  other 

0 K. 

coefficients  is  not  important.  A proof  of  this  result  is  given  in  Appendix  C. 
Let  us  multiply  both  sides  of  (22)  by  d [t+  (K-2i)Tl  and  integrate 

K T 1 

from  -®  to  * for  i=  0,1 K. 
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This  results  in  K + 1 simultaneous  equations,  which,  using  (10), 
become 


D 

D 

D 


K+1 

K+1 

K+1 


(0) 

(2t) 

(4t) 


D 


K+1 


D 


K+1 


D 


K+1 


(2t) 

(0) 

(2t) 


D (2K' 
K+r 


- 

- 

... 

1 

a 

1 

D^^,(2t)  ... 

K 

“2 

Dk^,(0|  ... 

• • 

• 

• 

• 

= 2n  t. 

I 

1 

• 

• 

* • 

• • 

^K+l<°)  _ 

^K 

Of 

K+1 

(23) 


where  we  have  substituted 
vironment  and 


a 


i+1 


tg  = 1 + r^  = 2 to  represent  the  marine  en- 

d (t+KT-2iT)dt,  i-0,l,2,...,K  (23a) 

K+ 1 


Substituting  for  D (e) 
Kt  1 


from  (21),  'f>e  find  that  (23)  reduces  to 


P(0) 

P(2t) 

P(4t).  . . P(2Kt) 

P(2t) 

P(0) 

P(2  r).  . . 

P(4t) 

P(2t) 

P(0)  ... 

R2K  ) 

• 

P(0) 

K , 

2 11  t 
1 ^ 


K+1 


(24) 


Note  that  the  (K  f 1 ) x (K  + 1)  matrix  on  the  left  has  the  Toeplitz 
structure.  The  terms  N(0),  N(2t),  . . . which  appear  in  P(0),  P(2t),  . . . , 
are  random  variables  with  known  statistics;  they  will  be  zero  if  the 
seismic  data  is  noise  free  (i,  e. , n(t)  -■  0).  Observe,  also,  that  the 


10 


first  and  last  elements  of  the  vector  on  the  left-hand  side  of  (24) 


are  unity  and  r , respectively,  by  virtue  of  the  property  which 
K 

we  just  stated  for  the  d,,  . causal  functional  equation. 

^ K+ 1 

Equation  (24)  provides  the  starting  point  for  our  procedure 

for  identifying  the  reflection  coefficients.  Note  that  in  general 

the  a are  functions  of  d (t),  a signal  wliich  is  not  determinable; 
i K+ 1 

liouever,  we  will  show  in  the  following  section  that  when  thi-  input  wave- 
let is  narrow  enough  (not  necessarily  a spik-e),  (24)  has  a unique 
solution  for  the  reflection  coefficients  in  terms  of  observable  data. 


7.  SPECIAL  CASE  OE  NARROW  WAVELET 


Let  us  now  consider  the  case  where  tlQ(t)  does  not  extend 


beyond  2t  , i.  e. 


dgft)  = 0 


t < 0 , t > 2t  • 


(25) 


Since  the  time  of  arrival  at  the  Kth  interface  is  Kt  and  the  time 
of  arrival  of  the  first  reflections  is  (K  + 2)t, 


t < Kt 


K 

= 2 n t.  d^ft -Kt)  KT<ts(K  + 2)T 

= more  complicated  terms  t>(K  + 2)T 


( 26a  ) 
( 26b  ) 
( 26c  ) 


If  this  condition  is  not  satisfied,  we  can  always  deconvolve  the 
data  to  achieve  this.  Since  the  requirement  here  is  not  to  deconvolve 
down  to  an  impulse  function  (only  (25)  has  to  be  satisfied),  this  re- 
sults in  a more  practical  solution. 
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From  (23a),  (26a),  (26b),  and  (25)  we  see  that 


+“> 

a.  = d (t)  d (t+ KT)dt  = 2 Tl  t.  ! d (t)  dt  , 


^ ■ J 


. 


K+1 


1 -'o  « 


and  that 


i+1 


d^(t)d^^^[t+Kr  - Zirldt  - 0 i 1 K 


(27) 


(28) 


Note  now  that  (24)  will  have  precisely  K + 1 unknowns,  K of  them  in 
the  vector  multiplying  the  Toeplitz  matrix  and  one  on  the  right-hand 
side,  OTj . 

Finally,  Normal  Equation  (24)  can  be  written  in  a compact 
matrix  form,  as 


K 


(2  9) 


where  P is  a (K+1)  X (K  + 1)  Toeplitz  matrix  with  the  first  row  being’^ 

K 

I.  P(0),  P(2t),  . . . , P{2Kt)];  a is  a K + 1 column  vector  with  first  and  last 

K. 

elements  1 and  r , respectively;  and,  C = col  (B,,,  0,  0,  . . . , 0)  and 
K.  K. 


K 


2 ^^2 
B = 2 n (1  - r ) r d (t)  dt 

^ ' -’o  ° 

The  Normal  Equation  (29)  can  be  solved  for  a . This  only  produces 

K 

one  of  the  K reflection  coefficients,  namely  r . We  will  show,  in  the 

K. 

following,  that  in  the  case  of  the  marine  environment,  nested  within  (29) 
are  a set  of  normal  equations,  the  solutions  of  which  produce  each  one 
of  the  reflection  coefficients.  The  absence  of  this  useful  property  in 
the  non-marine  case  renders  the  procedure  of  this  paper  inapplicable 


^'For  a narrow  wavelet  and  e > 2t,  the  calculation  of  P(e)  simplifies  since 
(18)  reduces  to 

_ 2 T 

P(e)  = -2  dt  . 
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in  that  case  [see  Appendix  Bl. 


I 

Let  us  now  hypothesize  a j -layer  system  (i.  e. , the  basement 
layer  is  the  i-th)  consisting  of  the  top  j-layers  of  the  above  K-layer 
system  (K^j).  Clearly,  from  (29),  we  have 

P.  a.  = C.  ( 30  ) i 

J J J 

where  a.  will  again  have  1 and  r^  as  first  and  last  elements.  We 

shall  now  show  that,  in  the  case  of  the  marine  environment,  is 

+ ^ (j  + 1 ) Toeplitz  matrix  composed  of  the  top  left  corner  of  J 

P : i.  e.  , its  first  row  is  given  by  [ P(0),  P(2t),  . . . , P(2jT)] . 

^ ! 
For  the  moment  let  us  ignore  the  additive  noise  term  in  (19). 

Let  us  denote  by  Uj(t)  the  response  of  the  j-layer  system  (i.  e. , the 

term  Uj(t)  in  Figure  2 is  replaced  by  u^j(t)).  In  (16),  due  to  the 

fact  that  dQ(t)  = 0 for  t>2T,  the  last  value  of  Uj(t)  contributing  to 

P(e)  is  Uj(2t  + C).  In  determination  of  P^  , with  elements  P(t), 

e = 0,  . . . , 2jT  , for  a j-layered  system,  therefore,  the  last  value  of 

_ ! 
u'j(t)  contributing  to  P^  is  u'^^  [ 2(j  + 1 )t].  On  the  other  hand,  Uj(t)  is 

the  response  of  the  K-layer  system,  and 

Uj(t)  = uj(t)  0 & t s 2(j+  1)T  ( 31  ) 

since  the  first  return  from  the  interfaces  below  the  j-th  will  not 
appear  earlier  than  t=2(j+l)T.  Hence,  the  elements  of  P^  , which 
are  functions  of  u'l  (t)  and  P which  are  functions  of  u,(t),  are 
identical  when  (31)  is  satisfied.  In  other  words,  the  numerical 
values  of  P(0),  . . . , P(2jT)  will  be  identical  to  those  of  the  K-layer 
system  for  all  j ^K.  Furthermore,  the  additive  noise  term  in  (19) 
is  independent  of  the  number  of  layers,  as  is  evident  from  (20). 


The  .-at  of  normal  equations  given  by  (30),  for  j = 1,,  . . ,K,  can 


now  be  solved  for  the  vectors  a.  and,  hence,  their  last  elements,  r , 

J j 

j = I , . . . , K.  The  matrix  P is  Toeplitz  and  consequently , the  Levinson 
algorithm  [11  can  be  used  to  solve  for  the  vectors  a.,  j = 1,..,,K 
recursively. 

Since  the  r.'s  are  reflection  coefficients,  for  the  solution  to  this 
J 

problem  to  be  acceptable,  each  r^  must  be  less  than  unity  in  magnitude. 

It  is  shown  in  Appendix  D that  any  solution  of  (30)  with  P.>0  for  all  j 

yields  a set  of  r.'s  which  satisfy  this  condition.  Moreover  if  P is 
J K 

positive  definite,  a compatible  solution  with  3.>0  is  guaranteed.  We 

see  therefore,  that  the  requirement  that  | r ^ | < 1 has  nothing  to  do  with 

a specific  method  of  solution  of  the  normal  equations  (i.  e.  the  Levinson 

procedure).  This  result  is  different  from  the  comparable  result  in  [1]- 

[3],  where  one  is  left  with  the  impression  that  a specific  method  of 

solution  leads  to  1 r . < 1 . 

' J 


[ 8.  REMARKS 

I 

1 

A comparison  between  the  procedure  of  this  paper  and  the 
standard  procedures  described  in  [11- [4]  is  warranted.  Let  us  consider 
the  noise  term  n(t)  to  be  white.  Clearly,  (20)  indicates  that  the  random 
variables  N(')have  finite  variances.  For  this  case  (n(t)  white),  had 
we  performed  the  necessary  deconvolution  and  sampling  required  by 
the  classical  approach  to  the  inverse  problem,  the  resulting  N(*) 
random  variables  would  have  infinite  variance,  clearly  rendering  the 

approach  meaningless.  Of  course,  "approximate"  deconvolution  will 

*If  n(t)  is  not  wide-band,  then  the  variance  of  N ( • ) may  not  be  infinite, 
but  will  be  very  large. 
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eliminate  this  problem  but  at  a great  sacrifice  in  the  information 


available  within  the  seismic  data.  It  should  also  be  noted  that,  for  the 
narrow  wavelets,  no  deconvolution  is  required  by  the  procedure 
outlined  in  this  paper. 

Finally,  the  procedure  of  this  paper  can  be  applied  (see  Appendix 
E)  to  the  classical  solution  of  the  inverse  problem  as  it  appears  in 

[i]-r4i. 

9.  CONCLUSIONS 

We  have  developed  a procedure  for  extracting  reflection 
(^Qgffidents  from  noisy  data  which  we  feel  is  a substantial  generalization 
of  similar  procedures  which  have  been  reported  in  the  literature. 
Associated  with  these  earlier  procedures  are  Standard  Assumptions 
and  Steps  (see  Introduction,  page  1)  which  include  requirements  that 
the  data  be  noise  free  and  that  the  observed  seismic  data  be  deconvolved. 
The  procedure  of  our  paper  avoids  these  restrictive  requirements. 
Furthermore,  our  procedure  totally  avoids  the  concepts  of  z-transforms, 
minimum  phase,  spectral  factorization,  forward  and  reverse  polynomial 
manipulations,  etc.  , which  appear  in  the  literature  on  this  subject. 
Finally,  since  our  derivation  is  so  straightforward,  it  suggests  a 
number  of  extensions,  including  the  following,  which  are  presently 
under  study:  (1)  nonstandard  locations  of  source  and  sensors  (e.  g.  , 
both  in  the  first  layer);  (2)  minimum  mean-square  estimation  of  the 
reflection  coefficients;  and  (3)  optimal  prefiltering  of  noisy  data 
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(exploiting  the  potential  utility  of  the  abstract  Kunetz  relation  in 
Appendix  A). 
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APPENDIX  A 


An  Abstract  Kunetz-Type  Relationship 

A generalized  Kunetz  relationship  was  derived  in  Section  4 
by  starting  with  the  following  state  equations, 


u.{t  + T)  = t.  u (t)  + r d (t) 

J J J J 

d^^l(.  + T)  = + t.d.(t) 


( A1  ) 


( A2  ) 


and  performing  the  following  operations:  (a)  time  shift  u.(t+T)  and 

d.  (t  + T)  by  e to  obtain  u.(t+T+e)  and  d (t  + T + e)  [this  is  done  by 
j+ 1 J 

operating  on  both  sides  of  (Al)  and  (A2)  by  a linear  advance  operator]; 

(b)  multiply  u.(t  + T)  by  u.(t  + T+E),  and  d (t+T)  by  d(t+T  + e);  (c)  add 
J J J J 

(t  /t  ) u (t  + T)  u, (t  + T + e)  to  d.(t+T)  d (t  + T+E);  (d)  integrate  the  resulting 

j j j J J J 

expression  from  -»  to  +»;  and  (e)  iterate  the  integrated  expression 
across  the  K-layers, 

Step  (c)  is  a pivotal  one;  for,  it  leads  to  a cancellation  of  terms 

common  to  both  (t./t!  ) u (t  + t)  u (t  + T + e)  and  d (t  + t)  d (t  + t + e).  Step 
3 J J J J J 

(e)  is  also  important  in  that  it  permits  us  to  use  the  boundary  condition, 
that 

Steps  (a),  (b),  and  (d)  involve  specific  operations.  We  show 

here  that  these  operations  can  be  abstracted  (i.  e.  , generalized),  and, 

that  we  can  obtain  an  abstract  Kunetz-type  relationship  which  is  valid 

for  much  more  general  operations  than  shift,  multiply,  and  integrate. 

Let  L denote  a linear  operator  for  the  j layer.  Operate  on 
3 

both  sides  of  (AI)  and  (A2)  with  L.  , to  show  that 
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L.  u^(t  + T)  = L.  + r.  L.  d.(t)  ( A3  ) 

L d (t  + T)  = -r.  L.  u.  (t)  + t.  L.  d.(t)  ( A4  ) 

j j+r  J J J+1  J J J 

Let  o denote  an  operation  that  satisfies  certain  group  properties  (e.  g. 
u could  be  multiplication,  convolution,  etc.  ) ; then, 

[L.  u.(t  + T)]  o u.(t+T)  = t'.^[L  u (t)]  o u (t) 

J j J J J J’’’ 

+ t.  r.  [ L.  u (t)]  o d (t) 

J J J J+1  J 

+ r^(L.  od.(t)  ( A5  ) 

J J J J 

and 

[L.dj^jd  + Tll  od.^,{t  + T)  = r^^[L.u.^j(t)l 

- 'i 

+ U\L.  d.(t))  o d.(t)  ( A6  ) 

J ^ J J J 

From  (A5)  and  (A6),  we  find  that 

[L.d  ,{t+T)]od.  ,(t  + T)  + (t./t'.)[  L u (t  + T)]  o u (t  + T) 

'•  J j+1  J + i J J J J J 

= f L.  u.^j(t)]o  u.^j(t)  + (t./L)[L.  d.(t)]  o d.(t)  ( A7  ) 

which  can  also  be  written,  as 

(L,d.^|(t,T))„d^^,(ttT)  - ou.^jdl 

. (t./t!  ) I [L,  d.(t)|  o d (t)  - [ L.  u (t  + T)]  o u (t+  T|) 
jj  Jj  J jj  J 


( A8  ) 


Let  N denote  another  linear  operator.  Operate  on  both  sides 


of  (A8)  with  N,  to  show  that 


= (t./t'.  )|Ni[  L.  d.(t)]  o d.(t)}  - n1[L.  u.(t  + T)]  ou.(t  + T)}  ] { A9  ) 

J J JJ  J JJ  J 

Now,  iterate  (A9)  backwards,  starting  with  j = K,  for  which  we 

know  that  that 


K 

n t. 

• 1 ’’ 

+ ° + = “k iNl[Lj  dj(t)]o  dj(t)} 

n t! 

. , 1 


- N|[Lj  Uj(t  + T)]o  Uj(t+T)}  } 

( AlO  ) 

Equation  (A  9)  is  an  abstract  Kunetz-types  relationship.  It 
reduces  to  Eq.  (12)  when 


and 


o = X 


all  j,  where  z is  an  e sec. 
advance  operator, 
(multiplication) 


00 

N = I ( • ) dt 

.CO 


( All  ) 


Equation  (AlO)  is  then  equivalent  to  Eq.  (13). 

The  applicability  of  the  abstract  Kunetz-type  relationship  using 
other  operations  [e.  g.  , = a low  pass  filter,  o - * (convolution), 

N = another  filtering  operation]  remains  to  be  studied.  What  is 


19 


interesting,  though,  is  the  fact  that  a very  general  relationship  - 


(A9)  or  (A  10)  - exists  for  a layered  media  system,  and,  that  this 
relationship  does  not,  in  general,  have  anything  to  do  with  energy 
or  an  energy  spectrum. 


APPENDIX  B 


Non-Marine  Case 

We  question  what  limits  the  results  of  this  paper  mainly  to  the 
marine  problem.  The  answer  is  two  fold.  First,  although  here  a 
structural  property  similar  to  (22)  is  still  valid  leading  to  (29),  the 
extension  of  (29)  to  (30)  is  only  possible  for  the  marine  case,  since,  in 
the  non -marine  case,  Eq.  (14)  is  d^(t+T)  = ‘ ^q^I 
that  the  cancellation  of  Uj(t)  u^(t+£)  terms  in  (15)  will  not  occur. 
Therefore,  P(e)  will  be  a function  of  the  entire  u^  (t)  for  tc  [O,®), 
regardless  of  how  narrow  dQ(t)  is.  This  suggests  that,  for  the  non- 
marine case,  it  seems  we  can  only  (at  least  directly)  solve  for  the  last 

reflection  coefficient  (r  ) and  not  r.,  j-1,...,  K-  1 . Second,  as  we  shall 

K J 

show  below,  the  nature  of  the  noise  term,  N(e),  is  more  complex  in  the 
non-marine  case,  to  the  point  that  its  statistics  cannot  be  obtained  from 
those  of  the  observation  noise  alone. 

Let 

PQ(e)  ^ 0^(6)  - ( B1  ) 

From  (10)  and  (11),  we  then  have 

^00  -f-OO 

PqIO  = J ^^(t)  dQ(t+  t)  dt  - J UQ(t)  UQ(t+  e)  dt  : ( B2  ) 

.CO  .00 

and,  from  (13), 


K 

n t. 


The  measurement  in  the  non-marine  case  is  given  by 


y(t)  = u (t)  f n(t)  ( B4  ) 


Since  ^^(t)  is  not  available,  we  may  use  y(t)  in  its  place  by 
defining  Pq(£),  as 


+<» 


i+CO 


^0^'^  = j ^^(t)  dQ(t  + e)  dt  - ^ ^ y{t)  y(t  + e)  dt  ( B5  ) 


which,  using  (B2),  can  be  written  as 


PqCe)  = PQ(e)  + N(e) 


( B6  ) 


where 


+®  +00 

N(e)  = - j n(t)  Up(t  + e)  dt  - J UQ(t)  n(t  + e)  dt 

-CD  .CD 

P+» 

- J n(t)  n(t  + e)  dt 


( B7  ) 


Finally,  (B3)  becomes 


K 

n t. 
n t! 


( B8  ) 


Note  that  here  the  statistics  of  N(e)  depend  on  the  statistics  of  n(t)  and 
the  availability  of  ^^(t):  however,  ^^(t)  is  not  available  and  can  only  be 
obtained  through  a noisy  measurement.  Additionally,  the  third  term  in 
N(e)  is  nonlinear  in  n{t).  Observe,  from  (20),  that  no  such  nonlinear 
noise  term  occurs  in  the  marine  case. 
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APPENDIX  C 


Proof  of  Structural  Property 

In  the  following,  a constructive  proof  of  Eq.  (22)  is  given. 
Equations  (1)  and  (2)  can  be  re-written  as 


'‘ji"  = l-l'i • 


( Cl  ) 


( C2  ) 


where  (Cl)  is  obtained  by  eliminating  the  term  dj(t)  between  (l)and 
(2).  For  j = 0,  (C2)  becomes 


do(t)  = 7'[*'o  + d (t+T)] 

0 


( C3  ) 


Replacing  for  Uj(t)  and  dj(t+T)  from  (Cl)  and  (C2)  yields 

^ tV  - t)  + rj  U2(t+  T)  + r^r^  d2(t)  + d2(t  + 2T)}  ( C4  ) 


Replacing  for  the  terms  in  the  right  of  (C4)  from  (Cl)  and  (C2)  yields 

= TTT:  1^0  + '’2  + 

0 1 2 

+ r^r^  d^(t  - T)  + (r^+r^r^)  d2(t+ T)  + d2(t  + 3t))  ( C5  ) 


Repeating  the  process  another  time  yields 

■*0"'  = TT^  + 

\J  I ^ J 

+ (^2'*' ^‘l'*' u^(t  + T)  + r^  u^(t+  3t)  + r^r^  d^(t  - 2t) 
+ ^1^2^  ^0*^2^  ^^4^*^  ^^l'^^l*'2^^2^3^  d^(t  + 2T) 


+ d^(t  + 4T) 


( C6  ) 
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Equation  (C3)  is  for  a one  layer  system  where  d^(t)  is  defined  at 
the  bottom  of  layer  1 (see  Figure  1),  That  equation  can  also  be  used  to 
describe  a zero  layer  system  (i.  e. , two  semi-infinite  half  spaces)  by 
setting  Uj(t)  = 0 and  moving  d^  from  the  bottom  of  layer  1 to  the  top  of 
that  layer.  This  is  accomplished  by  setting  t = t - t in  the  d^  term. 

The  resulting  equation  is 


todo(t)  = d^(t) 


(C7) 


We  proceed  similarly  for  (C4),  (C5),  and  {C6).  For  example,  (C6)  is 
for  a four  layer  system  where  d^(t)  is  defined  at  the  bottom  of  layer  4. 
That  equation  can  also  be  used  to  describe  a three  layer  system  by 
setting  u^(t)  = 0 and  moving  d^  from  the  bottom  of  layer  4 to  the  top  of 
that  layer.  This  is  accomplished  by  setting  t = t-T  in  all  d^  terms. 
The  resulting  equation  is 


n t.  dglt)  = d^(t+3T)  + ajd^(t+2T)  + a2d^(t-T)  + r^r^  d^(t-3T)  (C8) 


This  development  demonstrates  the  existence  of  our  equation  (22). 
Set  K=0  in  (22)  to  obtain  (C7),  and  set  K = 3 in  (22)  to  obtain  (C8). 


APPENDIX  D 


Some  Properties  of  Toeplitz  Matrices 


Consider  a set  of  a equations  of  the  type  encountered  in  Section  7: 


A.  X.  = b.  , 
1 


1 — Ifccv^n 


( D1  ) 


where  the  A^ , i-l,...,n  are  the  leading  principal  minors  of  the 

symmetric  Toeplitz  matrix  A = A , and  b.  and  x are  ixl  vectors 

n 1 I 

with  the  special  form 


with  b -K  , X = 1 (in  the  notation  of  Section  7,  A = P , x = a 
‘ ^ ^ i i i i 

^i“  ^i'  ^i"  ^i  ^i  ■ ^i^*  Several  properties  of  equations  of  the 
type  (Dl)  which  bear  on  the  solution  of  equation  (30)  will  now  be 


derived. 


First  note  that  if  A is  Toeplitz,  then  it  can  be  partitioned  as 


so  that  (Dl  ) is  equivalent  to  the  pair  of  equations 


A + m,  y.  = K. 


( D2  ) 


m.  + A y.  = 0 

i 1- 1 i 


( D3  ) 


Also,  if  det  A.  j ^ 0.  then  by  a well  known  determinant  identity  [8] 

' - 1 

det  A.  = det  A.  ,(A-m.  A.  ,m.) 

1 1-11  1 1- 1 1 

Furthermore,  if  det  A.  , ^ 0 ( D3  ) can  be  solved  for  y.  and  substituted 

i-l  1 

into  ( D2  ) yielding 

' -1 

K.  = A - m.  A.  , m. 

1 1 1 1-1  1 

Hence  if  det  A.  , ^ 0 we  have  the  identity 
1-  1 


det  A. 

K - ^ 

i det  A,  , 
1-1 


( D4  ) 


Theorem;  A symmetric  Toeplitz  matrix  A^  is  positive  definite  if 
and  only  if  for  each  i=l,...,n,  the  equation  (Dl)  has  a unique 
solution  y.  for  some  K.>0. 

Proof:  We  prove  sufficiency  by  induction  on  the  sequence  of  matrices 
A..  For  i=l  (Dl)  merely  implies  Suppose  now 

we  have  shown  that  Aj^>0,  i=l,...,k-l.  Applying  (D2)  for  i = k and 
the  fact  that  det  A > 0 by  the  induction  hypothesis  we  have  by  (D4) 

iC  • 1 

that 


det  A,  - K det  A,  , > 0 
k k k- 1 

which  completes  the  induction  step  since  Aj^>0  is  i^^  leading  principle 

minor  of  A,  . 

k 

To  prove  necessity  note  that  if  A>0  then  det  Aj^>0  for  i=l,...,n 
so  that,  for  each  i,  the  equations  (D2)  and  (D3)  have  a unique  solution 
with 
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> 0 


det  A. 


K 


det  A 


i-1 


Corollary:  If  A is  a symmetric  Toeplitz  matrix  with  leading  principal 

minors  A^ , i = 1 , . . . , n and  y[  = [ • • • 7^  ^1  ^ solution  of  (D1  ) for 


K >0,  then  1 y|  -^1  for  i - 1 n. 

1 1 


-1 


Proof:  By  our  theorem,  A.>0  so  that  Q.  = A.  >0.  Hence 

' 1 11 


[ 10  . . . 0 -1  ] Q. 


1 

0 


0 

-1 


Since  A^  is  symmetric,  Q.  is  symmetric  so  that 


Further, 


since  A^  is  Toeplitz 


■^1 

m 

■ Vi 

t ~ 

1 

m 

Vi. 

1' 

A 

1 -1 

so  that 


det  A 


i-1 


‘11 


det  A. 


= ‘*ii 


{ D5  ) 


Hence,  j j " ^ ® 


‘^il  < ‘111 


If  y^  satisfies  ( D1  ) withK^>0,  then  the  last  element  of  y^  satisfies 


y!"*  = q.,K.  < q,  ,K.  = 1 

’ll  1 11  i 
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since  by  (D4)  and  (D5),  1/K.. 

By  a similar  argument  we  can  show  that  qj  j f > 0 which 

implies 


- 28  - 


APPENDIX  E 


i 

I 

i 


i 


The  Classical  Discrete-Time  Solution  to  Inverse  Problem 

It  may  be  of  some  interest  to  relate  the  results  of  this  paper 
to  those  of  reference  [l]-[4]. 

Let  the  variable  t assume  discrete  values  which  are  a multiple 
of  T : i.  e. , 

t = kT  , k = 0, 1 , . ( El  ) 

Furthermore,  let  the  input  be  a unit  pulse  (unit  spike)  at  t = 0 

d^(t)  = A(0)  ( E2  ) 

With  this  notation,  the  only  change  to  be  made  in  the  derivation  is 
the  redefinition  of  quantities  Ej(e)  and  Uj(e)  in  (10)  and  (11).  Let 

e = iT  ( E3  ) 

A 

D.(kT)  = Z)  d.fkT]  d.[(k  + i)T]  ( E4  ) 

J i=-a,  J 3' 

A +® 

U.(kT)  = Z/  u.[kT]  u.[(k  + i)T]  ( E5  ) 

i=-oa  ■! 

The  energy  transfer  relation  (13)  remains  the  same,  and  (16)  is  re- 
placed by 

P(iT)  = 4 A(iT)  - 2 Uj(iT)  ( E6  ) 

We  also  disregard  the  presence  of  the  noise  term  in  (19).  Equation 
(22)  is  again  valid  for  *3^(1)  A^lt),  t = 0. 

The  remainder  of  the  derivation  follows  in  similar  manner  re- 
sulting in  the  following  normal  equation  (to  replace  (29))! 
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